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Abstract 


Exposure estimates inside space vehicles, surface habitats, and high altitude aircraft 
exposed to space radiation are highly influenced by secondary neutron production. The 
deterministic transport code HZETRN (High charge ( Z ) and Energy TRaNsport) has 
been identified as a reliable and efficient tool for such studies, but improvements to the 
underlying transport models and numerical methods are still necessary. In this paper, the 
forward-backward (FB) and directionally coupled forward-backward (DC) neutron 
transport models are derived, numerical methods for the FB model are reviewed, and a 
computationally efficient numerical solution is presented for the DC model. Both models 
are compared to the Monte Carlo codes HETC-HEDS (High Energy Transport Code - 
Human Exploration and Development of Space), FLUKA (FLUctuating KAskade), and 
MCNPX (Monte Carlo Neutral Particle extended), and the DC model is shown to agree 
closely with the Monte Carlo results. Finally, it is found in the development of either 
model that the decoupling of low energy neutrons from the light particle transport 
procedure adversely affects low energy light ion fluence spectra and the subsequent 
exposure quantities. A first order correction is presented to resolve the problem, and it is 
shown to be both accurate and efficient. 

1. Introduction 

Radiation exposure guidelines are a primary concern for the design of personal shielding, 
spacecraft, instrumentation, and missions. Consequently, there is significant interest in developing 
computational tools that allow shield analyses not only in simplified geometries, but in more complex final 
design models as well [Wilson et al. 2003]. The deterministic transport code HZETRN [Wilson et al. 1991, 
2006; Slaba et al. 2008], developed at NASA Langley Research Center has emerged in recent years as a 
reliable and efficient tool for such studies. It has shown reasonable accuracy in deep space Galactic Cosmic 
Ray (GCR), Solar Particle Event (SPE), and Low Earth Orbit (LEO) simulations when compared to either 
Monte Carlo results or experimental data [Wilson et al. 2005; Heinbockel et al. 2009]. However, such 
verification and validation has revealed a deficiency in the low energy neutron transport procedure [Shinn 
et al. 1994; Heinbockel et al. 2009]. HZETRN utilizes the straight ahead approximation in which all 
fragments are assumed to propagate in the same direction as the projectile. The assumption is accurate for 
high energy charged particles but begins to break down for low energy neutrons which are nearly isotropic 
[Alsmiller et al. 1965]. This is significant for heavily shielded space vehicles, surface habitats, and high 
altitude aircraft where secondary neutron production is important in exposure calculations [Getselev et al. 
2004], 

Several neutron transport models have been developed for HZETRN, with recent efforts focused 
on identifying an optimal bi-directional neutron transport model and solution method. The terms model and 
method are used extensively throughout this paper; model is used to refer to the set of governing transport 
equations, while method is used to refer to the analytic or numerical techniques used to solve a model. 
Heinbockel et al. [2000] and Clowdsley et al. [2000a, 2000b] developed the forward-backward (FB) 
neutron transport model. The FB model assumes that low energy neutrons can be split into forward and 
backward components, and multiple reflections from forward to backward (or vice versa) can be ignored. 
Feldman [2003] expanded on the work of Heinbockel and Clowdsley by developing the directionally 
coupled forward-backward (DC) neutron transport model. The DC model also assumes that low energy 
neutrons can be split into forward and backward components, but multiple reflections are accounted for in 
the governing transport equations. The numerical methods used for each model were significantly different 
as well. 

Heinbockel et al. [2000] and Clowdsley et al. [2000a, 2000b] used a multigroup method for 
solving the FB model. Multigroup methods have a long history in nuclear reactor theory and assume that 
fluences or cross sections are constant over small regions, or groups, of the energy spectrum [Marchuk and 
Lebedev 1986]. Feldman [2003] used a collocation technique along with first order finite differencing for 
the DC model. Collocation and finite differencing methods use polynomial expansions of the solution to 
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transform a differential equation into a system of linear equations [Deuflhard and Borne mann 2002]. In this 
case, the system was sufficiently large that it had to be solved numerically. 

Recently, Slaba et al. [2006, 2007] and Heinbockel et al. [2007] developed three methods for 
solving the FB model which they called the collocation method, the fixed-point series method, and Wilson's 
method. The latter two methods are both based on a Neumann series solution wherein each term of the 
series solves a simple set of equations related to the original model. This allowed for intensive verification 
of the multigroup method, and the collocation method was identified as the most accurate and 
computationally efficient of all the methods [Slaba 2007], Slaba [2007] also indicated that a combination of 
the methods could be used to obtain a computationally efficient Neumann series solution to the DC model. 
Finally, despite all of the attention given to neutron transport models and methods, accurately coupling 
these models back into HZETRN remained unresolved. The impact of such a coupling has not been 
previously studied in detail and must be examined if any of the neutron transport models are to be used in 
design studies. 

In this paper, we first give a summary of the neutron transport models and methods developed thus 
far and present an efficient method for solving the DC model. Next, comparisons are made between the FB 
and DC models and the Monte Carlo codes HETC-HEDS [Gabriel et al. 1995; Townsend et al. 2005], 
FLUKA [Fasso et al. 2003, 2005], and MCNPX [Briesmeister, 2000; MCNPX 2.6.0 Manual, 2008], 
Finally, the impact of decoupling low energy neutrons from the light ion transport procedure in HZETRN is 
also examined, and a first order coupling is presented. Fluence spectra and dose quantities are given to 
exhibit the accuracy of the proposed correction. 

2. HZETRN Description 

The one-dimensional Boltzmann transport equation for charged and neutral particles with the 
continuous slowing down and straight ahead approximations is given as [Wilson et al. 1991] 

B\4> j \ = E c jk {E,E')cl> k {x,E')dE' , (1) 


with the linear differential operator 


B M 


d_ 

dx 


— —— Sj (E) 
A j dE ] 




G ’ 


( 2 ) 


where (f>- is the fluence of type j ions at depth x with kinetic energy E (MeV/amu), A- is the atomic 
mass of a type /' particle, S-{E) is the stopping power of a type j ion with kinetic energy E , a j(E) is the 
total macroscopic cross section for a type j particle with kinetic energy E , and a- k {E, E ') is the 

macroscopic production cross section for interactions in which a type k particle with kinetic energy E ' 
produce a type j particle with kinetic energy E. Macroscopic cross sections are obtained by multiplying the 
corresponding microscopic cross section by the target particle mass density [Wilson et al. 1991]. Hereafter, 
it is assumed that all cross sections are macroscopic whether it is explicitly written or not. The summation 
over k on the right hand side of equation (1) will be explained shortly. 

Wilson et al. [1986, 1991, 2006] obtained approximate solutions to equation (1) by introducing 
the scaled quantities 


r ) = VjSpiE^faE), 

(3) 

s ]k (r,r') = S p (E)a jk (E,E'), 

(4) 


where S p (E) is the proton stopping power, r is the residual proton range 
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E dE 


( 5 ) 


r = f 

Jo Sp (E ') 


the scaling parameter v . = Z 2 / /T , and Z, is the charge of a type j ion. For neutrons, v is taken as unity 

in fluence scaling relations and zero in range scaling relations. This will be explained in more detail 
shortly. Equation (1) is now given in terms of the variables x and r as 



( 6 ) 


which can be inverted using the method of characteristics [Haberman 1998] and written as the Volterra 
type integral equation [Wilson et al. 2006] 

ipj (x, r) = 0,r- + v^x) 

E V „■ nx poo r ( r „,i\ (7) 

— / / e j{ 5 5 (r + v,x',r')^Ax - x\r') dr' dx' , 

, l/, J 0 J r+v-x J J 

k k J 

with the exponential attenuation term 




( 8 ) 


Note that v- in equation (3) and v- / v k in equation (6) are both fluence scaling relations, and v n = 1 in 
both cases to provide a nontrivial scaling. Conversely, wherever v appears as the argument of a fluence or 
cross section in equations (7) and (8), it is a range scaling relation and v = 0 . This convention is taken to 
reflect the absence of atomic interactions in neutron transport ( S n (E) = 0 ) and to avoid a trivial scaling for 
the neutron fluence. 

From here, the problem is split into two parts: heavy ions (A > 4) and light particles (.4 < 4) . 
For heavy ions, it is noted that projectile fragments have energy and direction very near that of the 
projectile (equal velocity assumption), while target fragments are produced nearly isotropically with low 
energy and travel only a short distance before being absorbed. This approximate decoupling of target and 
projectile fragments is discussed in detail by Wilson et al. [1991] and suggests that target fragments can be 
neglected in the heavy ion transport procedure (their contribution to dose is approximately accounted for 
after the transport procedure). The equal velocity assumption allows the production cross section to be 
recast as 


s j k ( r ’ r ') = v jk {r)8(r ~ r '), (9) 

or equivalently, in terms of energy as 

a jk(E,E') = a jk (E)S(E — E'). ( 10 ) 


Note that <Jj k (E) can be written interchangeably with E or r since, for a given E. equation (5) can be 

evaluated to obtain r. The absence of target fragments in the heavy ion transport procedure allows one to 
take the summation in equation (7) over all projectiles with mass greater than that of the fragment. If all 
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transported heavy ions are ordered according to mass, then the final marching procedure for heavy ions is 
given by [Wilson et al. 2006] 


ip Ax + h,r) = e ^ r ’ h ^ipAx,r + vh) 


k> j k 


jk 


h 

r + v.~ 
3 2 


A 


, ,h\ 

x,r + {v 1 +v k )-\ 


Ml / \h It. 

_|ft -a^r+ivj+v^-jh 


r + {v. + v,)- 



\ , h) 

3 \ 

\ 3 2) 


0(h 2 ) , 


(ID 


where h is the step size in g/cm 2 . The summation over all k > j reflects the fact that only heavy ion 
projectiles fragments are explicitly transported and that the transported heavy ions have been ordered 
according to mass. The upper summation limit in equation (11) can vary, but it is now common to use no 
fewer than 59 isotopes. See Cucinotta et al. [2006] for a discussion of isotopes grids. 

For light particles, both projectile and target fragments are included in the transport procedure; 
therefore, the summation is taken over all light particles. The broad energy distribution in collision events 
also indicates that equal velocity assumption cannot be used. The final marching procedure for light 
particles is given by [Wilson et al. 1991, 2006] 


ip Ax + h,r) = e ^ r ’ h ^ipAx,r + v-h) 


+ 


, „ V n 

E—f 

7A i 


l/ u J r+vJi/2 




x,r' 


v k 2 


dr' 


( 12 ) 


+0(h 2 


where the integrand has been simplified using 

F jk (r,r'-,h) = f Q s jk (r + u jZ , r ') dz , 


(13) 


and the neutron elastic interactions in this equation are handled separately as [Slaba et al. 2008] 


F^ l (r,r';h)=h{s* n (r,Ar')), 


(14) 


where F^fAr, r';h) denotes the elastic component and (s^ n (r,Ar')j is the average value of the 

differential neutron elastic cross section over a small interval, A r ' , centered at r'. 

The nature of the boundary condition, or environment, determines how equations (11) and (12) are 
evaluated and coupled. For most SPE environments, there is usually a negligible heavy ion component, and 
only equation (12) is evaluated. For GCR environments, equation (1 1) is evaluated for heavy ions; equation 
(12) is evaluated for light particles, and the summation term appearing in equation (11) is added. The 
summation is taken only over heavy ion projectiles, and it physically accounts for light ion production from 
heavy ion projectiles. The marching equations (11) and (12) will be referred to throughout this paper as the 
HZETRN marching algorithm, with the proper equations evaluated and coupled for a given environment. 

3. Neutron Transport Models 

To derive the neutron transport models, we first decompose the particle fluence and neutron 
production cross section into straight ahead and isotropic components. This approximation is justified by 
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Wilson et al. [1991] who showed that projectile fragments are produced with energy and direction very 
near that of the projectile, while target fragments are produced nearly isotropically with very low energy. 
For particle fluences, the decomposition is given as [Wilson et al. 2005] 

<f> j (x,E) = <f>f°(x,E) + <t>° a (x,E). (15) 


To decompose the neutron production cross section, we first distinguish between the nuclear 
reactive and elastic parts 


*nk(E,E') = a: k (E,E') + a* k (E,E') 


(16) 


where crfl (E, E ') = 0 for k ^ n . The angularly dependent nuclear reactive neutron production cross 
section can be represented as [ Wilson et al. 2005] 

a r nk (E,E',Sl,£l') = g{A T , E, 9) a' nk (E, E'), (17) 


where fl and Q' are unit vectors in the direction of the secondary neutron and projectile particle, 
respectively. The parameterization for g(A T ,E,9) is given by the piecewise definition [Ranft, 1980] 


g{A T ,E,9) 


Ne~ e2/X , 0 < 9 < tt/ 2 
^y e -7r 2 /4A ; tt/2 < 9 < 77 


(18) 


and A = (120 + 0.3671^,) / E , where E is the produced neutron kinetic energy in MeV, Aj, is the atomic 

mass of the target nucleus, 9 is the production angle with respect to the direction of propagation, and N is 
a normalization constant such that 


2t r J'J g(A r , E, 9) d9 = 1. 


(19) 


The parameterization constants, 120 and 0.36 have units MeV and MeV/amu, respectively, so that A is 
dimensionless. 

The nuclear reactive neutron production cross section can now be decomposed into forward 
(scattering angle 9 € [0, n/2 \ ) and backward (scattering angle 9 £ ( 77 / 2 , 77 ] ) components 


< k (E,E') = £)<,(£,£') + F(-)(^, £)<*(£,£') 

= ^ + \E,E') + a r Jr\E,E'), 

where the forward and backward coefficients are 


( 20 ) 


f (+) (^,£) = r /2 g(A T ,E,9)d9, (21) 

J 0 

(Aj, , E) = 1 - F {+ ) ( \ , E) . (22) 


The straight ahead and isotropic coefficients are therefore defined as [Clowdsley et al. 2000b] 
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F sa {A T ,E) = 1 — F lso (A T ,E), 

(23) 

F iso (A T ,E) = 2 FE)(A t ,E), 

(24) 


so that the nuclear reactive neutron production cross section can be expressed in terms of isotropic and 
straight ahead components 


<u(E,E') = F sa (Ap , E)cj r nk (E,E') + F iso (A, , E)a^ lk (E, E ' 


„r,sa 

® nk 


(E,E') + a^°(E,E'). 


(25) 


The decompositions in equations (20) and (25) are depicted in Figs. 1 and 2 for a 500 MeV proton 
projectile and an Aluminum target. It should be noted that in Fig. 2 the straight ahead component of the 
nuclear reactive neutron production cross section is dominant at high fragment energies, while the isotropic 
component is dominant at low fragment energies as discussed earlier. 

Equations (15), (16) and (25) are substituted into equation (1), and the HZETRN marching 
algorithm is used to solve the straight ahead transport equation 

B [ <f>? } = £ / ” 4 ° {E, E <) C {x,E')dE\ (26) 

k 

where the summation limits for equation (26) are environment dependent and have been discussed 
previously, and the production cross section K,- k (E, E ') is piecewise defined as 


*%(E,E') = 


vT(E,E') + a* k (E,E') ,3 = n 
a jk (E, E') , j ^ n 


(27) 



Figure 1 . Forward (+), backward (-), and total neutron nuclear reactive production cross section for a 500 

MeV proton projectile and Aluminum target. 
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Energy (MeV) 

Figure 2. Straight ahead (sa), isotropic (iso), and total neutron nuclear reactive production cross section for 

a 500 MeV proton projectile and Aluminum target. 



Figure 3. Neutron and 4 He fluence spectra at 50 g/cm 2 in an Aluminum target exposed to the February, 

1956 Webber SPE. 
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Note that for light particle projectile/fragment combinations, <j jk (E,E') will be fully energy dependent, 

while for heavy ion projectiles, a- k (E,E') incorporates the equal velocity assumption and assumes the 
simplified form in equation (10). 

The only difference between equations (1) and (26) is the form of the neutron production cross 
section; equation (1) uses the total nuclear reactive neutron production cross section, while equation (26) 
uses only the straight ahead component. Fig. 2 suggests that neutron fluences obtained from equation (26) 
will have attenuated low energy spectra, and an examination of the HZETRN marching algorithm indicates 
that a residual attenuation will occur in the low energy light ion spectra as well. This behavior is evident in 
Fig. 3. More will be said about this later. 

We are now left to solve for the isotropic neutron fluence 

B \€°\ = E fz a nk( E ,E')4 SO &E')dE' + n n {x,E), (28) 

k 

with the isotropic neutron source term 


V n (x,E) ee ^f”<T r Jr{ E , E ')W(x, E ') dE ' , (29) 

k 

The isotropic component of the charged particles is obtained by solving 

b [ <t>r ] = e 5e a * {E ’ E ') ^ so ^ E ') dE '> ( 3 °) 

k 

which will be handled later. For the moment, the summations in equation (28)-(30) are taken over all 
transported particles; the summations will be simplified later based on physical arguments. 

We now return to equation (28) and note that the contribution to the secondary neutron field from 
isotropic, low energy, charged target fragments will be small since the range of the range of these charged 
particles is small compared to their mean free path length (Wilson et al. 1991] as seen in Fig. 4. Thus, only 
the neutron projectile - neutron fragment integral need be retained. In this case, the isotropic neutron 
transport equation reduces to 

B[4>iT] = fz°nn( E , E ')<t>n°& E ')d E ' + V n (x, E )- OD 

The isotropic neutron fluence is now decomposed into forward and backward components 
according to 


<P:°(x,E) = cj>l(x,E) + . 


(32) 


The terms forward and backward have different meanings when applied to fluences and cross sections. The 
forward component of the isotropic neutron fluence refers to all neutrons propagating in the forward 
direction; the backward component of the isotropic neutron fluence refers to all neutrons propagating in the 
backward direction. The forward component of the neutron production cross section refers to those 
neutrons produced with scattering angle in the interval [0,7r/2] (no change in direction of propagation); the 
backward component of the neutron production cross section refers to those neutrons produced with 
scattering angle in the interval {i: 2, it] (change in direction of propagation from forward to backward or 
vice versa). 
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Figure 4. Range and mean free path length of 4 He in Aluminum. 


The nuclear reactive neutron production cross section is assumed to take the form of equation (20), 
thereby splitting it into forward and backward components. The elastic portion is also decomposed, but it is 
accomplished by considering the approximate relationship between scattering angle, pre-collision energy 
(. E ') and post-collision energy (E) given by Haffner [1967] 


E = E 


-j- 2A^, cos 9 1 

(4 + l) 2 


(33) 


Forward scattering occurs for 9 £ [0,7r/2] or for E' £ [E, E//3\, and backward scattering occurs for 
9 £ ( 77 / 2 , 7r] or for E' £ ( Ej(5 , E Ja\ where 


A t -if 

A + 1 , 

A\ + 1 

(A + !) 2 ' 


(34) 


(35) 


The neutron elastic production cross section is then decomposed into forward and backward components as 


T eZ(+) 


{E,E' 


< n (.E,E'), E'e[E,E/0] 
0 , otherwise 


(36) 


9 


a 


nn V ’ / 


a d n (E,E'), E' e(E/(3,E/a] 

0 , otherwise 


(37) 


The complete neutron production cross section can now be expressed as forward and backward components 
according to 

<r m ( E ’ E ') = + o^{E,E\ (38) 

with 

°%\ E , E ') = *T ] (E,E') + <(+)(£,£'), (39) 

*tt( E , E ') = + <«(£,£')■ (40) 


From Feldman [2003], substitution of equations (32) and (38) into equation (31) results in the DC neutron 
transport model 


B {+) W n ] = f” <tW (E,E') ft (x, E')dE' 
bH 1€} = f E *tt(E,E')<j> b n (x,E')dE' 


f E °° <tH (E, E ') ft (x, E')dE' + ri (x, E ) , 
(E,E') </>{ 0 v , E')dE'+r, b n (x,E), 


(41) 


(42) 


with the source terms r]l(x, E) and rj b (x,E) each taken to be one half of the isotropic source term 
r} (x,E) , and the linear differential operators are 


B M [<p\ = 

B ( ~\cj>} = 


7T + a " {E) 

OX 


- -”JE) 

OX 


<t>{x, E), 
<t>(x,E). 


(43) 

(44) 


Note that the minus sign in equation (44) accounts for a directional change in propagation. In equation (41), 
the second integral represents neutrons produced in the forward direction by backward propagating 
neutrons. In equation (42), the second integral represents neutrons produced in the backward direction by 
forward propagating neutrons. 

Slaba [2007] then showed that by making the simple approximation <f>^(x,E) m <p b (x,E) one 

obtains 


B (+) W n ] = /“ °nn (E, E ’) 4 (x, E')dE' + ri (x, E) , (45) 

B H [<t> b n ] = J7 v nn (E, E')cj> b n (x,E')dE' + r, b n (x, E) . (46) 

Equations (45) and (46) define the FB neutron transport model. The physical difference between the 
models can be deduced by examining the collision integrals. The DC model accounts for multiple changes 
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of direction (from forward to backward or vice versa) in the coupling source integrals, while the FB model 
only accounts for a single change of direction. That is, after the initial decomposition of the isotropic 
neutron source term into forward and backward components, it is assumed in the FB model that all forward 
neutrons propagate forward, and all backward neutrons propagate backward. The two models will be 
compared later in the paper. 

4. FB Neutron Transport Methods 

Heinbockel et al. [2000] and Clowdsley et al. [2000a, 2000b] used a multigroup method for 
solving the FB neutron transport model. The computational procedure is developed by partitioning the 
energy grid with respect to the elastic interaction parameter a defined in equation (34). Equations (45) and 
(46) are then integrated with respect to energy for some discrete set of energy values. The order of 
integration is switched in the collision integral and a mean value theorem is applied to evaluate the 
integrals. The result is an upper triangular system of first order ordinary differential equations solved using 
back substitution. Though the computational procedure turned out to be relatively efficient, the target 
dependent energy grid made energy grid convergence testing difficult, and selection of the free parameter 
in a mean value theorem came down to trial and error for several materials. 

Slaba et al. [2006, 2007] and Heinbockel et al. [2007] developed three different methods 
(collocation method, fixed point series method, and Wilson's method) for solving the FB model and found 
them to be in good agreement with the multigroup technique. The collocation method [Slaba 2007, 
Heinbockel et al. 2007] is implemented by assuming that forward and backward neutron fluences can be 
expressed as a sum of linear basis splines 


N 

3 = 1 
3=1 


( 47 ) 


(48) 


with 


B iW = 


(E - E^) / (E J - E^) , E £ [E-_ V E-] 

(e j+ 1 -e)/(e j+ 1 -e j ) , /•:e(/,>/T. l ; 

0 , otherwise 


(49) 


The approximations are substituted into the FB model resulting in two upper triangular systems of first 
order linear ordinary differential equations. Using back substitution, the analytic solution to the systems 
have been found to be [Slaba 2007, Heinbockel et al. 2007] 


^E n _ 3 ) = e a "-^ x) cj>l(0,E N _ .) + f’e'*->»-t*-’y n {x',E N _ J )dx' 
+(1 - k f 0 X e aN ^ N -^ X 'Ul(x\E N _ k )dx<, 


(50) 


P b n {x,E N _.) = e a »-^ {L - x) ^ n {L,E N _.) + f e^-^tiix^E^dx' 


+(! ~ 

k = 0 


j0J/_^ a N-j,N-k * e 


j 

x 


j{X '~ X U b n (x\E N _ k )dx\ 


(51) 
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where 6 is the Kronecker delta, L is the thickness of the material, and the coefficients a- are defined 

y 


% = + f™a,JE t ,E')B 3 (E')dE' . (52) 

The fixed-point series method [Slaba 2007, Heinbockel et al. 2007] is a Neumann series solution 
to the FB model. The forward and backward components are expressed as 

OO 

cj>l(x,E) = J2t\^E), (53) 

k = 1 

OO 

cj> b n (x,E) = J2bi k \x,E), (54) 

*=i 

where the k th zero order term satisfies the relations 

B (+) [jf ] = ZUx,E), (55) 

B { -% k) ] = e k - 1 (x,E), (56) 

and for k > 0 , the k th set of source terms are defined according to 

$ (*, E) = / ” <7 nn (77, 77 ')/o (fc) (*, 77 ') d77 \ (57) 

il (a-, 77) = / " (77, 77 (*, 77 ') dE ' . (58) 


The initial set of source terms ^ (x, E) and ^ (x, 77) are simply the forward and backward components of 
the isotropic neutron source term r] (x,E) . 

Wilson's method [Slaba et al. 2006, 2007; Heinbockel et al. 2007] is a combination of the 
collocation method and fixed-point series method. The forward and backward neutron fluences are 
expressed as the sum of a zero order and first order term by 

cj>f(x,E)= f 0 (x,E) + f^E), (59) 


<t>i(x,E) = b 0 (x,E) + \{x, E) . 


(60) 


Equations (59) and (60) are substituted into the FB model, and it is assumed that the zero order terms 
satisfy the relations 


B^[f 0 ] = V f n (x,E), 

(61) 

B(~%} = r, b n (x,E). 

(62) 
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The FB model is now reduced to 


5 (+) [/i] = J™cr nn (E,E')f l (x,E')dE' + t{(x,E), (63) 

bH 1\\ = f™a nn (E,E')b 1 (x,E')dE' + $(x,E) , (64) 


where the source terms (x,E) and ^(x. E) are calculated from the zero order terms by 


III 

ST 

fz°nn( E ,E%(x,E')dE', 

(65) 

£ l(x,E) = ^ 

f > OO 

f E <T m (E,E%(x,E')dE' . 

(66) 


Equations (63) and (64) are solved using the collocation method. 

Slaba et al. [2006, 2007] and Heinbockel et al. [2007] showed that the three methods produce 
almost identical results in a variety of shielding configurations exposed to both SPE and GCR 
environments and compared well with the multi group technique in all cases [Slaba 2007]. The work also 
revealed that the collocation method was the most computationally efficient of the four methods in terms of 
both memory requirements and run time. Hereafter, any neutron spectra labeled “FB model” are assumed to 
be generated using the collocation method. 

5. DC Neutron Transport Methods 

Feldman [2003] developed a numerical solution for the DC model. A collocation method was used 
to discretize the energy spectrum as in equations (47) and (48), and first order finite differencing was used 
to approximate the spatial derivatives and decouple the transport equations. Discretization of the energy 
and spatial domains in this way resulted in a large system of linear equations that was solved numerically 
using LU (Lower-Upper) factorization. LU factorization relies on writing the coefficient matrix as the 
product of a lower triangular and an upper triangular part so that forward and backward substitution can be 
used to obtain a solution [Demmel 1997]. It was reported that even single layer slab calculations quickly 
exhausted the memory capabilities of a desktop PC, and run times of over twenty hours were reported even 
on specialized computers with adequate memory [Feldman 2003]. 

We present here a computationally efficient Neumann series solution to the DC model based on 
the numerical methods verified for the FB model. The forward and backward components of the isotropic 
neutron fluence are expressed as a sum of zero order terms according to 


OO 

tfa,E) = £/o W (z,£), 


k = 1 


(67) 


OO 

tt(x,E) = '£bW(x,E) , 


k — l 


( 68 ) 


where the k th zero order term satisfies 

B W U f] = J^lV(E,E')f^\x,E')dE' + :,E), (69) 
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(70) 


B ( -¥o ] ] = / ^ fe) (ar, i? ^ ' + Cti {x,E), 

and can be solved using the collocation method outlined previously. Note that solutions to equations (69) 
and (70) include neutron production with scattering angle between 0 and 7^2 . Neutron production with 

scattering angle between 7^2 and 7r is accounted for in the k th set of source terms as 

i[{x,E) = f™a(J(E,E')bi k \x,E')dE', (71) 

e k (x, E) = f™ aH (E, E ')/ 0 (fc) {x,E')dE\ (72) 


for k > 0 . The initial set of source terms, ^ (x, E) and ^ ( x , E ) , are simply the forward and backward 

components of the isotropic neutron source term rj n (x, E) . 

While many terms in the Neumann series are required to adequately resolve the neutron fluence 
spectra (as many as 75 terms for each component), the computation is quite rapid since the coefficients a tj 

used in the collocation method and the cross sections appearing in the source terms in equations (7 1 ) and 
(72) can be pre-calculated and saved for repeated use. No published data exists from Feldman's work, and 
so any comparisons with data labeled “DC Model" were generated using the proposed Neumann series 
solution. 

6. Coupling to Light Ion Transport 

In the development of the neutron transport models, we did not solve the isotropic charged particle 
transport equation 


B [ \ = £ / “ a jk (E, E ') (*, E')dE' . (73) 

k 

At this point, the summation in equation (73) is taken over light particles only. This is justified by noting 
that is associated with low energy target fragments, and heavy ion target fragments are not transported 

in HZETRN due to their small range (Wilson et al. 1986, 1991, 2006]. 

The isotropic neutron fluence has already been obtained via one of the bi-directional neutron 
transport models, and so equation (73) can be expressed in the more explicit form 

B [ v;° \ = E f ” ** (E,E') or (x,E')dE' + Vj (x, E) , (74) 

k^n 

where the summation is taken over all charged, light particle projectiles {k ^ n ), and the neutron 
projectile term has been explicitly written as the source term 

/ OO 

E <r jn (E,E')<l>™{x,E')dE' . (75) 

We now approximate equation (74) by suppressing the remaining collision integrals. Such an 
approximation is expected to be valid since the mean free path of these low energy charged particles is 
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much larger than their range [Wilson et al. 1991 ] as was seen earlier in Fig. 4. In this case, equation (74) is 
now reduced to 




(76) 


for which an analytic solution exists and is easily found by using the method of characteristics [Haberman 
1998]. The solution was given by Wilson et al. [1975] in terms of the nuclear survival probability 


Pj{E) = exp 


-A 


[ E AAl dE . 1 

J o 


0 S/E< 


(77) 


and the residual energy E = R l [x + RAE)] as 


^°{x,E) = WW (f>?°(W 
3 Pj{E) S.{E) ^ K ^ 

rE AP(E') 

+ I 7 3 3 n [z + RAE) - RAE'),E'\dE\ 

Je P.{E)S\E) 3 3 jV 


(78) 


where Rj(E) is given by following the range energy relation 


w = a J! 


E dE' 


SAE ') 


(79) 


Equation (78) can be evaluated rapidly using numerical integration schemes and imposes small memory 
requirements. 


7. Results 


Solutions to equation (1) using the HZETRN marching algorithm defined by equations (11) and 
(12) will be labeled “SA” for straight ahead, and it is assumed that all particles (including neutrons) are 
transported in the forward direction. Solutions to equation (26) using the HZETRN marching algorithm 
with FB or DC neutrons added and no light ion coupling will be labeled “FB” or “DC”, respectively. 
Finally, solutions to equation (26) using the HZETRN marching algorithm with FB or DC neutrons and the 
suggested light ion coupling will be labeled “FBLI” or “DCLI”, respectively. The only difference between 
FB and FBLI or DC and DCLI is the presence of the light ion coupling term given in the previous section. 
Therefore, neutron fluences predicted by FB and FBLI (or DC and DCLI) will be identical. In all of the 
neutron fluence comparisons, FB and DC will be used, but could be interchanged with FBLI and DCLI 
without consequence. 

Results are now given at various depths in a 30 g/cm : water target behind a 20 g/cm 2 Aluminum 
shield exposed to the February 1956 Webber SPE spectrum [Quenby et al. 1959]. The Aluminum shield 
depth is indicative of spacecraft shielding, and the water depth is indicative of human tissue. This particular 
SPE spectrum was chosen because of the availability of Monte Carlo data. Fig. 5 shows the forward 
component of the neutron fluences at 20 g/cnr in the water target. The DC model compares quite well to 
HETC-HEDS and is within the bounds set by FLUKA and MCNPX; the FB and SA solutions predict much 
higher neutron fluences for energies below 100 MeV. Fig. 6 gives the backward neutron fluences at 20 
g/cnr 2 in the water target. The DC model again compares well to HETC-HEDS and is within the bounds set 
by FLUKA and MCNPX; it is a significant improvement over the FB model. There is no backward 
component of the SA model since all particles are transported in the forward direction. Finally, Fig. 7 gives 
the total neutron fluence at 20 g/cnr in the water target. 
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Figure 5. Forward neutron fluence at 20 g/cm 2 in a 30 g/cnT water target behind a 20 g/cm 2 Aluminum 
shield exposed to the February, 1956 Webber SPE. 



Figure 6. Backward neutron fluence at 20 g/cm 2 in a 30 g/cm 2 water target behind a 20 g/cm’ Aluminum 
shield exposed to the February, 1956 Webber SPE. 


16 



Figure 7. Total neutron fluence at 20 g/crrT in a 30 g/cnr water target behind a 20 g/cm’ Aluminum shield 

exposed to the February, 1956 Webber SPE. 


As in each of the previous figures, the DC model compares well to HETC-HEDS and is within the 
bounds set by FLUKA and MCNPX. Most importantly, it now appears that FIZETRN compares as well to 
the Monte Carlo codes as they compare with each other. Figs. 8-10 give the proton, 3 H, and 4 He fluence 
spectra generated from the SA, FB, DC, FBLI, and DCLI models at 20 g/cm 2 in the water target. First, note 
that the spectra predicted by the FB and DC models are identical. This behavior is expected since neither of 
these models include the suggested light ion coupling, and the charged particle fluences are determined by 
solutions to equation (26). 



Figure 8. Proton fluence at 20 g/cm 2 in a 30 g/cm 2 water target behind a 20 g/cm’ Aluminum shield 

exposed to the February, 1956 Webber SPE. 
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Figure 9. 3 H fluence at 20 g/cnr in a 30 g/cm 2 water target behind a 20 g/cm 2 Aluminum shield exposed to 

the February, 1956 Webber SPE. 



Figure 10. 4 FIe fluence at 20 g/cm 2 in a 30 g/cm’ water target behind a 20 g/cm 2 Aluminum shield exposed 

to the February, 1956 Webber SPE. 


Next, the spectra predicted by the FBLI and DCLI models are similar but not identical. In fact, the 
FBLI model seems to over-estimate the DCLI model in all three cases. This trend is caused by the light ion 
coupling source term in equation (74) and the neutron transport model used. If the FB model is used for 
neutron transport, then the light ion source term in equation (74) will be over predicted and so will the 


18 


isotropic light ion fluences. Surprisingly, the FBLI and DCLI light ion spectra agree quite well even though 
the FB and DC neutron fluences in Fig. 7 are almost an order of magnitude different at the lowest energies. 

Most importantly, Fig. 9 and Fig. 10 show that a significant reduction in low energy fluence 
spectra can occur if the suggested light ion coupling term is not included. Though little difference is seen in 
the proton fluence spectra (Fig. 8), the 4 FIe spectrum drops by almost a factor of four if the light ion 
coupling term is not included. Similar results are also seen for the other light ions transported in HZETRN 
( 2 FI and 3 He). Such deviations can alter integrated quantities as well. For example, Tables 1 and 2 give the 
dose and dose equivalent values at various depths in the water target. As expected, these values are all 
reduced if the light ion coupling is not included. Dose and dose equivalent values were not available from 
FLUKA and MCNPX at the time of this report. 


Table 1. Dose (cGy) in a 30 g/cm 2 water target behind 20 g/cm 2 of Aluminum exposed to the February, 
1956 Webber SPE spectrum. 


Depth 

g/cm 2 

SA 

DCLI 

FBLI 

DC 

FB 

HETC- 

HEDS 

0 

6.97 

6.99 

6.98 

6.95 

6.93 

7.20 

5 

3.78 

3.78 

3.78 

3.75 

3.75 

4.04 

10 

2.28 

2.28 

2.29 

2.26 

2.26 

2.50 

20 

0.98 

0.98 

0.99 

0.96 

0.97 

1.10 

30 

0.49 

0.48 

0.49 

0.47 

0.47 

0.53 


Table 2. Dose equivalent (cSv) in a 30 g/cm 2 water target behind 20 g/cm 2 of Aluminum exposed to the 
February, 1956 Webber SPE spectrum. 


Depth 

g/cm 2 

SA 

DCLI 

FBLI 

DC 

FB 

HETC- 

HEDS 

0 

11.70 

11.83 

11.50 

11.36 

11.03 

8.87 

5 

6.21 

6.18 

6.27 

5.80 

5.84 

4.85 

10 

3.94 

3.87 

4.01 

3.56 

3.69 

3.03 

20 

1.87 

1.79 

1.94 

1.59 

1.71 

1.36 

30 

1.05 

0.91 

1.07 

0.78 

0.91 

0.63 


8. Conclusions 

A directionally coupled neutron transport model and computationally efficient solution method 
were presented. Results compared well with HETC-HEDS and were within the bounds set by FLUKA and 
MCNPX. Therefore, the neutron transport algorithm in HZETRN now agrees with the Monte Carlo codes 
to the extent that they agree with each other. A first order approximation was also presented for coupling 
the neutron models to low energy light ion transport, and results indicated that such a coupling is indeed 
necessary not only for accurate estimates of low energy fluence spectra, but for integrated quantities such 
as dose and dose equivalent as well. It is recommended that future studies with HZETRN and bi-directional 
neutrons use the directionally coupled neutron transport model along with the first order light ion coupling 
presented here. 
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